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Abstract 

We investigate the liquid-solid transition of two dimensional hard 
spheres in the presence of gravity. We determine the transition tem- 
perature and the fraction of particles in the solid regime as a function 
of temperature via Even-Driven molecular dynamics simulations and 
compare them with the theoretical predictions. We then examine the 
configurational statistics of a vibrating bed from the view point of the 
liquid-solid transition by explicitly determining the transition tem- 
perature and the effective temperature, T, of the bed, and present a 
relation between T and the vibration strength. 

P.A.C.S: 05.20-y, 64.70Dv, 05.20.Dd, 51.10+y. 



The hard sphere model has been quite successful in explaining the macro- 
scopic properties of dense fluids, or gases from the microscopic point of view 
[1]. At the molecular level, the potential energy of the hard spheres due to 
gravity is small in comparison to the thermal fluctuations and it has usually 
been ignored. However, when the mass of the constituent particle is macro- 
scopic in quantity, as in the case of granular materials [2], gravity cannot be 
ignored. The purpose of this Letter is to demonstrate the existence of a grav- 
ity induced liquid-solid phase transition of hard spheres. This transition is 
an intrinsic transition associated with any system where the excluded volume 
interaction is dominant. Such a system cannot be compressed indefinitely, 
and must exhibit a coherent low energy state. In the hard sphere system, 
gravity introduces a potential energy, and each available site is associated 
with an energy state. Then, the formation of a solid at the bottom below 
the transition point is nothing but a massive occupation of the low energy 
state at the low temperature, which is the Fermi gas in metals, the Bose con- 
densate in the two dimensional quantum Hall systems [3] , the energy storage 
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mechanism into a single state for biological systems [4], and a mechanism 
to produce coherent light in the context of lasers [5], and the liquid-solid 
transition in a hard sphere system under gravity, which is the subject of 
the current work. We will determine via Even-Driven Molecular Dynamics 
simulation the transition point and the thickness of the boundary layers as 
a function of external parameters, and make a careful comparison with the 
theoretical predictions [6]. Next, a new and nontrivial by-product of our in- 
vestigation is to view the configurational statistics of the vibrating bed [7,8] 
from the view point of the liquid-solid transition of hard spheres. This will 
certainly help one to compare the configurational statistics and other ther- 
modynamic properties of vibrating beds with those equilibrium properties of 
hard sphere systems. 

Transition temperature and the thickness of boundary layers: Consider a 
collection of elastic hard spheres of mass m and diameter D, confined in a two 
dimensional (x, z) container with an open boundary at the top. Gravity acts 
along the negative z direction. The system is in contact with the thermal 
reservoir at a temperature T in such a way that the average kinetic energy of 
each hard sphere, T = m < v 2 +v 2 > /2 with < .. > being the configurational 
average. We now start from T = 0, at which point all the particles are 
essentially in a solid regime, and the density profile is simply a rectangle. If we 
gradually increase the temperature, fluidization starts from the surface, and 
the boundary layers appear. One may estimate the thickness of the boundary 
layer, h, by a simple energy balance between the kinetic and potential energy: 
mgDh ~ \m < v 2 >~ T. From this, one may obtain the size of the solid-like 
regime, or equivalently its dimensionless height, say (f(T): 

where \x is number of layers of the original rectangle. Eq.(l) predicts the 
existence of a critical temperature, T c , at which point a phase transition from 
a one phase(solid) to a two phase regime (liquid-solid) occurs. By setting, 
(f(T c ) = 0, we find the mean field result: T C M F = fimgD. Since the boundary 
layer exists only when both phases coexist, T c must be the temperature at 
which point the system becomes fully fluidized. One may equally define the 
critical temperature as a point at which the density at the bottom layer, <f> , 
becomes the closed packed density <f> c , i.e: (f> (T c ) = C . We now rewrite 
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eq. (1) in terms of the critical temperature, and recast the size of the solid 
region, in term of T/T c , as 



frCOAx = (i 



T 

c 
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A more precise estimate of the transition temperature was given in ref. [6] 
within the framework of Enskog theory [9]. In particular, the following ex- 
pression for the density profile, 0(C), was obtained as a function of the di- 
mensionless variable z/D: 

—f3(( — p) = ln(f> + Ci4> + c 2 /og(l — a<f>) + c 3 / (1 — a<f)) + c 4 /(l — a(p) 2 (3a) 

with the constant (5p given by: 

(3p = lri(p + ci<p + c 2 Zn(l - a<p ) + c 3 / (1 - a ) + c 4 /(l - a<p ) 2 = f(<p ) (36) 

where f3 = mgD/T, and c\ = 2a 2 /a 2 f ~ 0.0855, c 2 = —^(a 1 — 2a 2 /a)/a 2 m 
—0.710 c 3 = — c 2 ,c 4 = — «i/a + Q; 2 /a 2 /)« « 1.278. Note that the relation 
between the volume fraction v and <j) is given by: v = tt(D/2) 2 N/V = 7T0/4. 
If one integrates the density profile, and imposes the sum rule, Jq°° 4>(Qd( = 
4> ^i then one finds that this sum rule breaks down at the temperature T c , 



The departure from the mean field theory is the appearance of a factor <f> / /i 
in (4), where O = (4/7r)(7r/2 v / 3) = 1.154700538.. and fi = 111.52274... Eq. 
(2) remains unchanged. We now present MD data to test Eq.(2) and (4). 

Molecular Dynamics simulations of gravity induced liquid-solid transition: 
We have used the Event Driven(ED) Molecular Dynamics code, and refer the 
readers to references [10] for details of the algorithm regarding the collision 
dynamics that take into account the rotation of hard spheres, and a way 
to handle the inelastic collapse. The thermal reservoir of our system was 
modeled using white noise driving [11], which kicks each particle so that 
the average kinetic temperature of each particle is the same as that of the 
reservoir, and hence, the kinetic temperature of the system. Note that we 
are not driving the system by connecting the bottom wall to the temperature 
reservoir, which was often used as a model for a vibrating bed. 



T c = mgDfKpo/ fj, 



(4) 
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We present in Fig.l a typical configuration below the transition temper- 
ature (T < T c ), at which about 17 layers of particles condense and form a 
crystal near the bottom (Fig. la). More precisely, the particles first form a 
loose hexagonal crystal and progressively evolve into a compact hexagonal 
lattice structure. The solid line in the density profile(Fig.lb) is the Enskog 
profile given by (3a), which was shifted to fit the data beyond the crysal 
regime. We point out here that (i) this shift is not an arbitrary parameter, 
but should be uniquely chosen to fit the data, (ii) this shift in fact determines 
the measured size of the solid by simulations. The density in the solid regime 
is then fit by a straight line as shown in the figure. The oscillations in the 
solid regime are real, but it is simply the finite size effect, i.e, the hexagonal 
packing in a finite lattice has two more particles in alternative layers. This 
oscillation must disappear in the thermodynamic limit. 

The critical temperature T c is determined as the temperature at which 
point a compact hexagonal crystal is formed from the bottom layer, beyond 
which point, the density at the bottom layer remains constant at <p = 1.15, 
and this hexagonal structure is permanelty retained. We point out that a 
loosely hexagonal crystal forms at a temperature, T c ', which is somewhat 
larger than T c . Between T c and T c ', particles squeeze themselves, expelling 
holes, and progressively forming a compact hexagonal crystal. Note that a 
few vacancies created during this crystallization do not anneal but stay in the 
system (Fig. la). Now, in order to carry out the quantitative analysis of the 
formation of a crystal beyond the transition temperature, we have measured 
the size of the solid as mentioned above, namely by shifting the Enskog profile 
(i.e. Fig. lb), and plotted it at different temperature T < T c as a function 
of the scaled variable T/T c for 1000 particles of m = 2.090 • 10~ 6 , D = 
0.001m and /i = 20. The solid line in Fig.2a is the prediction Eq.(2). The 
excellent agreement between the theory and simulations is a confirmation of 
(i) the existence of the gravity induced liquid-solid transition of hard spheres, 
and (ii) the validity of the suggested mechanism of this transition via the 
disappearance of particles from the liquid and their settlement into the solid 
regime as predicted by Enskog theory [6]. 

Next, we present our new analysis of the vibrating bed from the view point 
of the liquid-solid transition discussed above. It has been fairly well estab- 
lished that the configurational statistics of the vibrating beds seem identical 
to the equilibrium statistics of a molecular gas at an equal packing fraction [8], 
yet the relation between the vibrational strength, T, and the corresponding 
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equilibrium kinetic temperature has remained largely undetermined. There 
has been a previous attempt to relate T to the Fermi temperature [12], which 
is not the same as the kinetic temperature, but essentially the compactiv- 
ity [13]. In the present work, we will establish a specific relation between 
the vibration strength and the kinetic temperature, and test its validity via 
simulations. 

At a low vibration strength, experimental data [7] seems to clearly in- 
dicate two distincitve regimes: solid regime near the bottom where there 
are very little particle movements, and the liquid regime near the surface 
where particles are dynamically active exchanging their positions via col- 
lisions. Hence, the system presented in re.f [7] is below the liquid -solid 
transition temperature. We will determine both the transition tempera- 
ture, T c , and the effective temperature of the system, and then measure the 
size of the solid region and compare it with the prediction made by Eq.(2). 
The control parameters are given in ref. [7], namely the particle diameter 
D = 2.99mm, and the dimensionless initial layer thicnkness /i = 10.2, from 
which we determine the normalized critical temerature of the vibrating bed 
T c /mg = (xD(f) / fi Q = 0.607mm. The effective temerature of the system is 
then determined by fitting the tail region and shifting the Enskog profile, by 
( . We find, T/mg = 0.36mm, and ( Q = 4.41 layers from which we mea- 
sure the size of the solid as z Q = <f) ( /D = 12mm m 4.0 layers, while the 
predicted dimensionless height of the solid region, £>, is: (f = [J>(l — T/T c ) m 
4.15 layers. The previous fitting of the density profile by the Fermi 

profile was also satisfactory, but was found to be most difficult near the 
rounded region, which the Enskog profile fits quite well. One advantage of 
the present method of analyzing the configurational statistics of the vibrat- 
ing bed might be that the global kinetic temperature can now be associated 
with the vibrating bed, and hence comparison can be made between the ex- 
perimentally determined configurational statistics of the vibrating bed and 
those of the hard spheres in thermal contact with the heat reservoir. The 
specific relation between the two can be obtained by comparing the thermal 
expansion of the hard spheres and the kinetic expansion of the vibrating bed. 
The thermal expansion is simply the increase in the center of mass Az(T), 
which can be computed by the Enskog profile near the tail, and the solid 
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rectangle. We find: 



Az(T) 



D/i ( 2\Ai\<f>, 
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where the constant |Ai| = O Jo[f(p4>o) — f(4>o)] [p4>of'(p4>o)dp] = 5503.531806 
with f(x) given in (3a). Note that the correction is second order in T. Let 
H (T) be the single ball jump height on the surface [16]. Then, by equating 
Az(T) and H (T)g/uj 2 , we find the desired relation: 



where u is the vibration frequency. Putting all the values, eq.(6) predicts 
T/T c = 0.663, which is close to the measured value of 0.593 above. 

In conclusion, two points are in order. First, we have demonstrated in 
this Letter that the point at which the Enskog description of hard spheres 
fails indeed signals the liquid-solid transition, and such a failure arises via 
the breakdown in the particle conservation. The missing particles form a 
condensate at the bottom, which essentially determine the fraction of parti- 
cles in the solid regime, and in turn the thickness of boundary layers. Since 
only a fraction of grains are mobilized under shear [14], and avalanches and 
many interesting dynamics occur in these thin boundary layers [2], such a 
determination should be of technological importance. Second, since Enskog 
theory is a truncation of BBGKY [15] hierarchy at the third order, the ex- 
istence of gravity induced liquid-solid transitions of hard spheres must have 
some interesting consequences to higher order kinetic theory, in particular 
with regard to the dynamic behaviors. Unlike particles in the liquid regime, 
those particles in the solid regime are largely confined in cages and fluctuate 
around fixed positions. Their motions resemble the lattice vibrations rather 
than binary collisions, and it may be a little peculiar, albeite not unphysical, 
to attempt to describe the lattice vibrations by the kinetic theory. If so, such 
a description must include much more than binary collisions. Hence, it is 
not unphysical to see that these particles disappear from the kinetic equa- 
tion at the level of the Enskog approximation. However, as discussed in the 
beginning and demonstrated in this paper, this gravity induced liquid-solid 
transition is not a peculiar phenomenon associated with Enskog equation, but 
rather an intrinsic transition inherent in a system where an excluded volume 
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interaction is dominant. The formation of a solid at the bottom is the ap- 
pearance of a massive occupied low energy state due to the Pauli exclusion 
principle. Therefore, the breakdown in the sum rule, the necessary shift of 
the density profile due to the formation, and its upward spread of the closed 
packed regime should persist because the Pauli exclusion principle is in action 
in real space, even if one may use different approximations [16-18] or may try 
a different form for the pressure, such as the form suggested by Percus-Yevick 
[19], and/or in higher order truncation. It only disappears in the limit when 
the closed volume packing density, v becomes one, which is possible only in 
the case of an ideal Appolonian packing [20]. Finally, we point out that the 
presence of dissipation does not alter the condensation picture at all [21], 
if the velocity distribution remains Gaussian. Recent experiments [22] have 
demonstrated the non-Gaussian nature of the velocity distribution, but if the 
dissipation is small, which is the case for the simulations carried out in this 
work, the deviation from Gaussian should be small. 
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Figure Captions 

Fig. 1. (a) Snap shot at T < T c , where about 17 layers form a crystal, (b) 
The fitting of the density profile is the combination of the Enskog profile 
(Eq.3) and the rectangle ( straight line). 

Fig. 2. The fraction of the hard spheres in the condensed regime as a function 
of T/T c with N = 1000, \i = 20, g = 981cm/sec 2 , and m = 1.047 • 10" 6 . 
(square). The data points are obtained by uniquely determining the shifting 
position of the Enskog profile, and the solid line is the prediction Eq.(2). 
Fig. 3. Experimental density profile of the granular materials in a vibrating 
bed (ref.5). The fit was by the Enskog profile near the surface, and the 
rectangle below £>. 
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